
#include "itkOrientedImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkAddImageFilter.h"
#include "itkShiftScaleImageFilter.h"
#include "itkThresholdImageFilter.h"
#include "itkExtractImageFilter.h"
#include "itkResampleImageFilter.h"
#include "itkLinearInterpolateImageFunction.h" 
#include "itkCenteredAffineTransform.h" 
#include "itkFlipImageFilter.h"
#include "itkIdentityTransform.h"
#include "itkRecursiveGaussianImageFilter.h"
#include <iostream>
#include <fstream>

#define PI (3.141592653589793)


	typedef float    PixelType;
	const unsigned int      Dimension = 3;
	typedef itk::OrientedImage< PixelType, Dimension >         ImageType;
	typedef itk::ImageFileReader< ImageType >        ReaderType;
	typedef itk::ShiftScaleImageFilter< ImageType, ImageType > ShiftScaleFilterType; 
	typedef itk::AddImageFilter<ImageType, ImageType, ImageType > AddFilterType; 
	typedef itk::ThresholdImageFilter< ImageType > ThresholdFilterType; 
	typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractFilterType; 
	typedef itk::FlipImageFilter< ImageType > FlipType;

	// PET/CT: Comparison of Quantitative Tracer Uptake Between Germanium and CT Transmission Attenuation-Corrected Images 
	// http://jnm.snmjournals.org/cgi/content/figsonly/43/9/1137
	const float bilinear_scale_hu[] = {-1000,0,3071};
	const float bilinear_scale_mu[] = {0,0.093,0.262};

	// Compute the y intercept of each of the bilinear scales
	// y = m * (x + h )
	// mu = m * ( hu + h )
	// h =  mu / m - hu
	
	// Compute constants for the ShiftScaleImageFilter

	const float m1 = (bilinear_scale_mu[1]-bilinear_scale_mu[0])/(bilinear_scale_hu[1]-bilinear_scale_hu[0]); // The slope
	const float shiftScale1[] = {
	bilinear_scale_mu[1]/m1-bilinear_scale_hu[1], // The y intercept	
	m1
	};
	const float m2 = (bilinear_scale_mu[2]-bilinear_scale_mu[1])/(bilinear_scale_hu[2]-bilinear_scale_hu[1]); // The slope
	const float shiftScale2[] = {
	bilinear_scale_mu[1]/m2-bilinear_scale_hu[1], // The y intercept
	m2
	};

	const float ADD_FILTER_DELTA = 1e-8; // This is used when adding bone and tissue images together
